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The  NASA  STI  Program  Office  ...  in  Profile 


Since  its  founding,  NASA  has  been  dedicated  to 
the  advancement  of  aeronautics  and  space 
science.  The  NASA  Scientific  and  Technical 
Information  (STT)  Program  Office  plays  a  key  part 
in  helping  NASA  maintain  this  important  role. 

The  NASA  STI  Program  Office  is  operated  by 
Langley  Research  Center,  the  Lead  Center  for 
NASA's  scientific  and  technical  information.  The 
NASA  STI  Program  Office  provides  access  to  the 
NASA  STI  Database,  the  largest  collection  of 
aeronautical  and  space  science  STI  in  the  world. 
The  Program  Office  is  also  NASA's  institutional 
mechanism  for  disseminating  the  results  of  its 
research  and  development  activities.  These  results 
are  published  by  NASA  in  the  NASA  STI  Report 
Series,  which  includes  the  following  report  tj^es: 

•  TECHNICAL  PUBLICATION.  Reports  of 
completed  research  or  a  major  significant 
phase  of  research  that  present  the  results  of 
NASA  programs  and  include  extensive  data 
or  theoretical  analysis.  Includes  compilations 
of  significant  scientific  and  technical  data  and 
information  deemed  to  be  of  continuing 
reference  value.  NASA's  counterpart  of  peer- 
reviewed  formal  professional  papers  but 
has  less  stringent  limitations  on  manuscript 
length  and  extent  of  graphic  presentations. 

•  TECHNICAL  MEMORANDUM.  Scientific 
and  technical  findings  that  are  preliminary  or 
of  specialized  interest,  e.g.,  quick  release 
reports,  working  papers,  and  bibliographies 
that  contain  minimal  annotation.  Does  not 
contain  extensive  analysis. 

•  CONTRACTOR  REPORT.  Scientific  and 
technical  findings  by  NASA-sponsored 
contractors  and  grantees. 


•  CONFERENCE  PUBLICATION.  Collected 
papers  from  scientific  and  technical 
conferences,  symposia,  seminars,  or  other 
meetings  sponsored  or  cosponsored  by 
NASA. 

•  SPECIAL  PUBLICATION.  Scientific, 
technical,  or  historical  information  from 
NASA  programs,  projects,  and  missions, 
often  concerned  with  subjects  having 
substantial  public  interest. 

•  TECHNICAL  TRANSLATION.  English- 
language  translations  of  foreign  scientific 
and  technical  material  pertinent  to  NASA's 
mission. 

Specialized  services  that  complement  the  STI 
Program  Office's  diverse  offerings  include 
creating  custom  thesauri,  building  customized 
data  bases,  organizing  and  publishing  research 
results  . . .  even  providing  videos. 

For  more  information  about  the  NASA  STI 
Program  Office,  see  the  following: 

•  Access  the  NASA  STI  Program  Home  Page 
at  http://www.sti.nasa.gov 

•  E-mail  your  question  via  the  Internet  to 
help@sti.nasa.gov 

•  Fax  your  question  to  the  NASA  Access 
Help  Desk  at  301-621-0134 

•  Telephone  the  NASA  Access  Help  Desk  at 
301-621-0390 

•  Write  to: 

NASA  Access  Help  Desk 

NASA  Center  for  AeroSpace  Information 

7121  Standard  Drive 

Hanover,  MD  21076 
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Summary 

The  Weibull  distribution  has  been  widely  adopted  for  the  statistical  description  and  inference  of 
fatigue  data.  This  document  provides  user  instructions,  examples,  and  verification  for  software  to  analyze 
gear  fatigue  test  data.  The  software  was  developed  presuming  the  data  are  adequately  modeled  using  a 
two-parameter  Weibull  distribution.  The  calculations  are  based  on  likelihood  methods,  and  the  approach 
taken  is  valid  for  data  that  include  type  1  censoring.  The  software  was  verified  by  reproducing  results  pub¬ 
lished  by  others. 


Introduction 

This  document  provides  user  instructions,  examples,  and  verification  data  for  software  to  analyze  gear 
fatigue  test  data.  The  software  was  developed  on  the  presumption  that  the  data  are  adequately  modeled 
using  a  two-parameter  Weibull  distribution.  It  is  based  on  the  likelihood  methods  described  by  Meeker 
and  Escobar  (1998).  The  software  can  be  used  to  determine 

1 .  Maximum  likelihood  estimates  of  the  Weibull  distribution 

2.  Data  for  contour  plots  of  relative  likelihood  for  two  parameters 

3.  Data  for  contour  plots  of  joint  confidence  regions 

4.  Data  for  the  profile  likelihood  of  the  Weibull  distribution  parameters 

5.  Data  for  the  profile  likelihood  of  any  percentile  of  the  distribution 

6.  Likelihood-based  confidence  intervals  for  parameters  and/or  percentiles  of  the  distribution 

The  software  can  account  for  tests  that  are  suspended  without  failure.  The  statistical  terminology 
for  suspended  tests  is  “censoring.”  The  analytical  approach  for  the  software  is  valid  for  type  I  censoring, 
which  is  the  removal  of  unfailed  units  at  a  prespecified  time.  Confidence  regions  and  intervals  are  calcu¬ 
lated  using  the  likelihood  ratio  method.  Guidance  for  the  philosophy  and  interpretation  of  statistical  confi¬ 
dence  intervals  can  be  found  in  the  text  of  Hahn  and  Meeker  (1991). 

The  Weibull  distribution  has  been  widely  adopted  for  the  statistical  description  and  inference  of 
fatigue  data.  The  Weibull  cumulative  distribution  function  F  is  defined  as 


F{t)  =  1  -  exp 


where  t  is  the  time  to  failure,  P  is  the  shape  parameter,  and  T|  is  the  scale  parameter.  For  the  Weibull 
distribution,  the  probability  density  function  can  have  a  variety  of  shapes,  and  the  hazard  function  may 
be  an  increasing  or  decreasing  function  of  time.  The  shapes  of  these  curves  depend  only  on  the  shape 
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TABLE  I.— DATA  SET  1 
[Data  from  Krantz  et  al.,  2000.] 


Test 

number 

Test  time, 
millions  cycles 

Test  result 

1 

40.26 

Fail 

lire 

2 

53.94 

3 

59.88 

4 

67.68 

5 

95.76 

6 

134.22 

7 

198.48 

8 

256.20 

1 

9 

299.52 

Suspended 

10 

301.56 

'  11 

303.66 

12 

304.08 

13 

305.80 

14 

306.90 

15 

335.40 

1 

TABLE  IL— DATA  SET  2 


[Data  from  Townsend  and  Shimski,  1994.] 


Test 

number 

Test  time, 
millions  cycles 

Test  result 

1 

7.08 

Failure 

2 

8.15 

3 

15.52 

4 

21.34 

5 

33.37 

6 

35.70 

7 

38.31 

8 

40.26 

9 

40.74 

10 

49.47 

11 

51.70 

12 

55.29 

13 

59.17 

14 

70.18 

15 

127.07 

16 

130.95 

17 

133.38 

r 

18 

300.00 

Suspended 

19 

300.00 

Suspended 

20 

300.00 

Suspended 

parameter.  Typically,  surface  fatigue  data  for  gears  have  distributions  with  shape  parameters  less  than 
three.  Therefore,  the  probability  distribution  functions  typically  are  skewed  right.  For  a  given  Weibull 
distribution,  one  can  determine  the  mode,  mean,  median,  variance,  and  higher  moments  analytically  by 
employing  the  gamma  function  (Cohen,  1973). 

The  software  developed  herein  was  written  in  the  FORTRAN  90  programming  language.  Users  of  the 
software  are  required  to  write  short  main  programs  that  provide  the  data  to  be  analyzed  and  also  call  the 
appropiiate  subroutines.  The  examples  provided  in  this  manual  were  compiled  using  the  Microsoft 
PowerStation  Fortran  Compiler,  version  4.0  (Professional  Edition). 

The  data  analyzed  as  examples  for  this  document  (tables  I  and  II)  are  those  from  gear  fatigue  experi¬ 
ments.  These  data  have  previously  been  analyzed  using  regression-based  methods,  and  the  results 
reported  (Krantz  et  al.,  2000).  The  data  are  again  analyzed  and  reported  in  this  document  using  software 
based  on  likelihood  methods. 

Verification  information  for  the  software  is  provided  in  appendix  A.  The  software  was  verified  by 
reproducing  results  published  by  Meeker  and  Escobar  (1998).  Appendix  B  provides  descriptions  and 
instructions  for  using  subroutines.  Some  calculations  require  the  use  of  routines  from  a  standard  math 
subroutine  library  (IMSL,  1997;  IMSL  is  a  registered  trademark  of  Visual  Numerics,  Inc.)  or  an  equiva¬ 
lent.  Subroutine  wmaxll  (appendix  B)  is  based  on  the  code  of  Keats,  Lawrence,  and  Wang  (1997). 


Estimates  of  Distribution  Parameters 

Data  Set  1 

The  main  program  that  analyzed  the  data  of  table  I  to  provide  maximum  likelihood  estimates  of  the 
Weibull  distribution  shape  and  scale  parameters  is  presented  below  as  source  code  listing  1.  Allocatable 
dimensioning  is  used  (lines  8  and  9)  to  make  certain  that  the  array  sizes  match  the  value  of  “n,”  the  vari¬ 
able  describing  the  number  of  tests  (line  7).  For  the  array  “censor,”  a  value  of  “zero”  is  given  to  indicate  a 
test  run  to  failure  whereas  a  value  of  “one”  is  given  to  indicate  a  test  suspended  without  failure  (lines  10 
to  15).  A  call  to  subroutine  wmaxll  provides  maximum  likelihood  estimates  of  shape  and  scale  parameters 
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(line  40).  With  the  maximum  likelihood  values  for  the  shape  and  scale  available,  subroutine  wblinverse 
provides  maximum  likelihood  estimates  for  any  desired  percentile  of  interest  (lines  48  to  52). 

The  main  program  produced  the  output  presented  below  following  source  code  listing  1 .  The  maxi¬ 
mum  likelihood  estimates  for  shape  and  scale  parameters  are  1.037  and  376.7,  respectively. 


Source  Code  Listing  1 

Line  Source  Line  Microsoft  Fortran  PowerStation  Compiler.  Version  4.0 

1  1  main  program  to  analyze  superfinish  gear  fatigue  test  results 

2  implicit  none 

3  integer  n, i 

4  real,  allocatable : :  time ( ; ) 

5  integer,  allocatable:;  censor (: ) 

6  real  shape, scale, frac, life 

7  n  =  15 

8  allocate (time (n) ) 

9  allocate (censor (n) ) 

10  do  i=l,8 

11  censor (i)  =  0 

12  end  do 

13  do  i-9,n 

14  censor (i)  =  1 

15  end  do 

16  open (20,  f ile=" superdat2 . txt" ) 

17  !  HERE  WE  NEED  TO  ENTER  THE  ARRAY  OF  TEST  TIMES 


18 

time  (1)  = 

40.26 

19 

time (2 ) = 

53.94 

20 

time (3 ) = 

59.88 

21 

time (4) - 

67.68 

22 

time (5) = 

95.76 

23 

time (6) = 

134.22 

24 

time (7) = 

198.48 

25 

t ime ( 8 ) = 

256.20 

26 

time (9) = 

299.52 

27 

time (10) = 

301.56 

28 

time (11) = 

303.66 

29 

time ( 12 ) = 

304.08 

30 

time (13 ) = 

305.80 

31 

time (14 ) = 

306.90 

32 

time (15)= 

335.40 

33 

do  i=l,n 

34  if  (censor (i)  .eq.  0  )  then 

35  write (20, 100)  time(i) 

36  else  if  (censor (i)  .eq.  1  )  then 

37  write (20 , 101)  time(i) 

38  end  if 

39  end  do 

40  call  wmaxll (time , censor, n, shape , scale) 

41  write (20, 105) 

42  write (20, 106)  shape,  scale 

43  100  format (f 10 .2, '  test  failed  ') 

44  101  format (f 10 .2, '  test  suspended  ') 

45  105  f ormat ( '  maximum  likelihood  estimates  ') 

46  106  f ormat ( '  shape  =  ',el2.4,'  scale  =  ',el2.4) 

47  frac  =0.0 

48  do  i=l,9 

49  frac  =  frac  +  0.1 

50  call  wblinverse (shape, scale, frac, life) 

51  write (20 , 107)  frac, life 

52  end  do 

53  107  format (f 6 . 2 , '  proportion  life  =  ',fl0.4) 

54  close (20) 

55  stop 

56  end 
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Output  produced  by  this  main  program  was  as  follows: 


40.26  test  failed 
53.94  test  failed 
59.88  test  failed 
67.68  test  failed 
95.76  test  failed 
134.22  test  failed 
198.48  test  failed 
256.20  test  failed 
299.52  test  suspended 
301.56  test  suspended 
303.66  test  suspended 
304.08  test  suspended 
305.80  test  suspended 
306.90  test  suspended 
335.40  test  suspended 
maximum  likelihood  estimates 
shape  =  .1037E+01  scale  =  .3767E+03 

.10  proportion  life  =  43.0405 

.20  proportion  life  =  88.7278 

.30  proportion  life  =  139.4505 

.40  proportion  life  =  197,1552 

.50  proportion  life  =  264.6016 

.60  proportion  life  =  346.2899 

.70  proportion  life  =  450.5648 

.80  proportion  life  =  596.0463 

.90  proportion  life  =  841.8331 


Estimates  of  Distribution  Parameters 

Data  Set  2 


The  main  program  that  analyzed  the  data  of  table  II  to  provide  maximum  likelihood  estimates  of  the 
Weibull  distribution  shape  and  scale  parameters  is  presented  below  as  source  code  listing  2.  The  source 
code  listing  2  minors  source  code  listing  1 ,  differing  only  in  the  lines  defining  the  number  and  results  of 
experiments  (lines  7  to  32).  The  main  program  produced  the  output  presented  below  following  source 
code  listing  2.  The  maximum  likelihood  estimates  for  shape  and  scale  parameters  are  0.8616  and  103.0, 
respectively. 


Source  Code  Listing 


Source  Line  Microsoft  Fortran  PowerStation  Compiler.  Version  4.0 


1  !  main  program  to  analyze  superfinish  gear  fatigue  test  results 

2  implicit  none 

3  integer  n, i 

4  real,  allocatable : :  time ( : ) 

5  integer,  allocatable:;  censor{:) 

6  real  shape, scale, frac, life 

7  n  =  20 

8  allocate (time (n) ) 

9  allocate (censor (n) ) 

10  do  i=l,17 

11  censor (i)  =  0 

12  end  do 

13  do  1=18, n 

14  censor (i)  =  1 

15  end  do 

16  open (20,  f ile=" superdat 1 . txt" ) 

17  I  HERE  WE  NEED  TO  ENTER  THE  ARRAY  OF  TEST  TIMES 

18  time(l)=  7.08 

19  time(2)=  8.15 
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20  time(3)=  15.52 

21  time(4}=  21.34 

22  time(5)=  33.37 

23  time(6)=  35.70 

24  time(7)=  38.32 

25  time(8)-  40.26 

26  time(9)=  40.74 

27  time(10)=  49.47 

28  time{ll)=  51.70 

29  time(12)=  55.29 

30  time(13)=  59.17 

31  time(14)=  70.81 

32  time (15)=  127.07 

33  time(16)=  130.95 

34  time(17)=  133.38 

35  time(18)=  300.00 

36  time(19)=  300.00 

37  time(20)=  300.00 

38  do  i=l,n 

34  if  (censor (i)  .eq.  0  )  then 

39  write (20, 100)  time(i) 

40  else  if  (censor(i)  .eq.  1  )  then 

41  write  (20 , 101)  tiine(i) 

42  end  if 

43  end  do 

44  call  wmaxll (time , censor, n, shape , scale) 

45  write (20, 105) 

46  write (20, 106)  shape,  scale 

47  100  format (f 10.2, '  test  failed  M 

48  101  format (f 10 . 2 , '  test  suspended  ') 

49  105  f ormat ( '  maximum  likelihood  estimates  ') 

50  106  f ormat ( '  shape  =  ',el2.4,'  scale  =  ',el2.4) 

51  frac  =  0.0 

52  do  i=l, 9 

53  frac  =  frac  +  0.1 

54  call  wblinverse (shape, scale, frac, life) 

55  write (20, 107)  frac, life 

56  end  do 

57  107  format (f6 .2 , '  proportion  life  =  ',fl0.4) 

58  close (20) 

59  stop 

60  end 

Output  produced  by  this  main  program  was  as  follows: 

7.08  test  failed 
8.15  test  failed 
15.52  test  failed 
21.34  test  failed 
33.37  test  failed 

35.70  test  failed 
38.31  test  failed 
40.26  test  failed 
40.74  test  failed 
49.47  test  failed 

51.70  test  failed 
55.29  test  failed 
59.17  test  failed 
70.81  test  failed 

127,07  test  failed 
130.95  test  failed 
133.38  test  failed 
300.00  test  suspended 
300.00  test  suspended 
300.00  test  suspended 
maximum  likelihood  estimates 
shape  =  .8616E+00  scale  =  ,1030E+03 

.10  proportion  life  =  7.5571 

.20  proportion  life  =  18.0554 

.30  proportion  life  =  31.1180 
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.40 

.50 

.60 

.70 

.80 

.90 


proportion 

life  = 

47.2135 

proportion 

life  - 

67.2835 

proportion 

life  = 

93.0217 

proportion 

life  = 

127.7066 

proportion 

life  = 

178.8617 

proportion 

life  = 

271.0445 

Joint  Confidence  Regions 

Data  Set  1 

The  main  program  that  analyzed  the  data  of  table  I  to  determine  joint  confidence  regions  for  the 
simultaneous  estimates  of  the  shape  parameter  and  10-percent  life  is  presented  below  as  source  code  list¬ 
ing  3.  Output  files  were  written  to  provide  data  for  likelihood  ratios  and  confidence  levels  (lines  25  to  28). 
Allocatable  dimensioning  is  used  (lines  30  and  31)  to  make  certain  that  the  array  sizes  match  the  value  of 
“n,”  the  variable  describing  the  number  of  tests  (line  29).  For  the  array  “censor,”  a  value  of  “zero”  is 
given  to  indicate  a  test  run  to  failure  whereas  a  value  of  “one”  is  given  to  indicate  a  test  suspended  with¬ 
out  failure  (lines  32  to  37).  The  range  of  values  for  the  shape  and  scale  parameters  must  be  provided  (lines 
59  to  62).  A  call  to  subroutine  wjlone  provides  the  data  required  for  contour  plots.  The  subroutine  writes 
results  to  FORTRAN  units  numbered  9  and  10.  The  output  from  the  program  is  two  data  files,  each  con¬ 
taining  1681  lines.  Each  line  consists  of  three  numbers  corresponding  to  a  10-percent  life  value,  a  shape 
parameter  value,  and  either  the  likelihood  ratio  or  the  likelihood-ratio-based  joint  confidence  value.  The 
joint  confidence  values  are  displayed  as  a  contour  plot  (fig.  1),  providing  a  graphical  display  of  plausible 
values  for  simultaneous  estimates  of  the  10-percent  life  and  the  shape  parameter. 


Figure  1 . — Likelihood-ratio-based  joint  confidence  regions  for  data  of  table  I. 
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Line  Source  Line  Microsoft  Fortran  PowerStation  Compiler.  Version  4.0 

1  I  a  main  program  to  calculate  likelihood  based  probability  contours 

2  1  and  relative  likelihood  contours  of  Weibull  distributions  from  censored  data 

3  1 

4  !  this  program  presumes  that  the  two  parameters  to 

5  !  be  plotted  are  the  shape  parameter  and  the  10%  life 

6  !  version  1.0 

7  !  written  by  Tim  Krantz  On  2-11-2001 

8  implicit  none  !  declare  all  variables 

9  integer  n  1  the  number  of  tests 

10  real,  allocatable : :  time(:)  1  the  test  times 

11  integer,  allocatable::  censor{:)  1  the  censoring  information 

12  real  tenmax, tenmin  I  the  max  and  min  values  for  the  10  percent  life 

13  !  to  be  used  in  calculations 

14  real  shapemin, shapemax  i  the  max  and  min  values  for  the  shape  factor 

15  1  to  be  used  in  calculations 

16  integer  i  I  a  do  loop  counter 

17  !  change  program  as  needed  for  a  particular  case 

18  !  change  only  lines  between  the  two  lines  of  stars 

19  I  *************************************************************** 

20  1  the  outputs  for  this  program  will  be  written  to  these  two  files 

21  1  the  file  formats  are  as  follows  (each  has  3  columns) 

22  I  column  1  is  the  shape  parameter 

23  1  column  2  is  the  10  percent  life 

24  i  column  3  is  the  likelihood  ratio  or  the  confidence  number  as  a  percentage 

25  1  file  for  unit  #10  is  the  confidence  number  data 

26  1  file  for  unit  #9  is  the  likelihood  ratio  data 

27  OPEN  (10,  FILE  =  "super-prob.txt") 

28  OPEN  (9,  FILE  =  "super-ratio.txt") 

29  n  =  15 

30  allocate (time (n) ) 

31  allocate (censor (n) ) 

32  do  i=l,8 

33  censor (i)  =  0 

34  end  do 

35  do  i=9,n 

36  censor (i)  =  1 

37  end  do 

38  1  HERE  WE  NEED  TO  ENTER  THE  ARRAY  OF  TEST  TIMES 

39  time(l)=  40.2600 

40  time (2) =  53.9400 

41  time(3)=  59.8800 

42  time(4)=  67.6800 

43  time(5)=:  95.7600 

44  time(6)=  134.220 

45  time(7)=  198.480 

46  time(8)=  256.200 

47  time(9)=  299.520 

48  time(10)=  301,560 

49  time(ll)=  303.660 

50  time(12)=  304.080 

51  time(13)=  305.800 

52  time(14)=  306.900 

53  time(15)=  335.400 

54  I  here  we  set  the  bounds  for  10  percent  life  and  shape  factors 

55  !  tenmax  is  the  largest  10  percent  life  value 

56  I  tenmin  is  the  smallest  10  percent  life  value 

57  I  shapemax  is  the  largest  shape  value 

58  I  shapemin  is  the  smallest  shape  value 

59  tenmax  =  140. 

60  tenmin  =  2.0 

6 1  shapemax  =  2 . 

62  shapemin  =  0.4 

63  !  none  of  the  lines  below  need  to  be  changed  to  run  a  particular  case 

64  1  ******** ********************** ****************** ****** ************** 

65  1  the  call  to  this  subroutine  checks  for  errors  in  the  censoring  array 

66  call  ccheck (n, censor) 

67  I  now  we  ask  for  the  solution 

68  call  wj lone (time, censor, n, shapemax, shapemin, tenmax, tenmin) 
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Data  Set  2 


69  close (10) 

70  close(9) 

71  stop 

72  end 


The  main  program  that  analyzed  the  data  of  table  II  to  determine  joint  confidence  regions  for  the 
simultaneous  estimates  of  the  shape  parameter  and  10-percent  life  is  presented  below  as  source  code  list¬ 
ing  4.  The  source  code  listing  4  min*ors  source  code  listing  3.  The  output  from  the  progi'am  is  two  data 
files,  each  containing  1681  lines.  Each  line  consists  of  three  numbers  con*esponding  to  a  10-percent  life 
value,  a  shape  parameter  value,  and  either  the  likelihood  ratio  or  the  likelihood-ratio-based  joint  confi¬ 
dence  value.  The  joint  confidence  values  are  displayed  as  a  contour  plot  (fig.  2),  providing  plausible 
values  for  simultaneous  estimates  of  the  10-percent  life  and  the  shape  parameter. 


Figure  2. — Likelihood-ratio-based  joint  confidence  regions  for  data  of  table  il. 


Source  Code  Listing  4 

Line  Source  Line  Microsoft  Fortran  PowerStation  Compiler.  Version  4.0 

1  I  a  main  program  to  calculate  likelihood  based  probability  contours 

2  i  and  relative  likelihood  contours  of  Weibull  distributions  from  censored  data 

3  I 

4  !  this  program  presumes  that  the  two  parameters  to 

5  !  be  plotted  are  the  shape  parameter  and  the  10%  life 
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6  !  version  1.0 

7  !  written  by  Tim  Krantz  On  2-11-2001 

8  implicit  none  !  declare  all  variables 

9  integer  n  !  the  number  of  tests 

10  real,  allocatable : :  time ( : )  1  the  test  times 

11  integer,  allocatable::  censor (: )  !  the  censoring  information 

12  real  tenmax, tenmin  !  the  max  and  min  values  for  the  10  percent  life 

13  !  to  be  used  in  calculations 

14  real  shapemin, shapemax  1  the  max  and  min  values  for  the  shape  factor 

15  !  to  be  used  in  calculations 

16  integer  i  !  a  do  loop  counter 

17  I  change  program  as  needed  for  a  particular  case 

18  I  change  only  lines  between  the  two  lines  of  stars 

19  1  *************************************************************** 

20  !  the  outputs  for  this  program  will  be  written  to  these  two  files 

21  J  the  file  formats  are  as  follows  (each  has  3  columns) 

22  !  column  1  is  the  shape  parameter 

23  !  column  2  is  the  10  percent  life 

24  I  column  3  is  the  likelihood  ratio  or  the  confidence  number  as  a  percentage 

25  I  file  for  unit  #10  is  the  confidence  number  data 

26  1  file  for  unit  #9  is  the  likelihood  ratio  data 

27  OPEN  (10,  FILE  =  "base-prob.txt") 

28  OPEN  (9,  FILE  =  "base-ratio.txt") 

29  n  =  20 

30  allocate (time (n) ) 

31  allocate (censor (n) ) 

32  do  i^l,17 

33  censor (i)  -  0 

34  end  do 

35  do  i=18,n 

36  censor (i)  =  1 

37  end  do 

38  !  HERE  WE  NEED  TO  ENTER  THE  ARRAY  OF  TEST  TIMES 

39  time(l)=  7.08 

40  time(2)=  8.15 

41  time(3)=  15.52 

42  time(4)=  21.34 

43  time(5)=  33.37 

44  time (6) =  35.70 

45  time (7) =  38.32 

46  time(8):=  40.26 

47  time(9)=  40.74 

48  time (10)=  49.47 

49  time(ll)=  51.70 

50  time(12)=  55.29 

51  time(13)=  59.17 

52  time(14)=  70.81 

53  time(15)=  127.07 

54  time (16)=  130.95 

55  time(17)=  133.38 

56  time(18)=  300.00 

57  time(19)=  300.00 

58  time(20)=  300.00 

59  !  here  we  set  the  bounds  for  10  percent  life  and  shape  factors 

60  !  tenmax  is  the  largest  10  percent  life  value 

61  !  tenmin  is  the  smallest  10  percent  life  value 

62  !  shapemax  is  the  largest  shape  value 

63  !  shapemin  is  the  smallest  shape  value 

64  tenmax  =  30. 

65  tenmin  =  1.5 

6  6  shapemax  =  1.4 

67  shapemin  =  0.4 

68  !  none  of  the  lines  below  need  to  be  changed  to  run  a  particular  case 

g9  t  ******************************************************************** 

70  1  the  call  to  this  subroutine  checks  for  errors  in  the  censoring  array 

71  call  ccheck (n, censor) 

72  1  now  we  ask  for  the  solution 

73  call  wj lone (time, censor, n, shapemax, shapemin, tenmax, tenmin) 

74  close  (10) 

75  close (9) 

76  stop 

77  end 
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Profile  Likelihood  and  Profile-Likelihood-Based  Confidence  Intervals 


Data  Set  1 


The  main  program  that  determined  the  profile  likelihood  and  profile-likelihood-based  confidence 
intervals  for  the  10-percent  life  for  the  data  of  table  I  is  presented  below  as  source  code  listing  5.  The 
main  program  deteimines  the  profile  likelihood  over  a  range  of  values  for  a  specified  percentile  of  inter¬ 
est,  in  this  case  the  10th  percentile.  Any  percentile  can  be  investigated  by  changing  one  line  in  the  code 
for  variable  “pfrac”  (line  49).  Allocatable  dimensioning  is  used  (lines  18  and  19)  to  make  certain  that  the 
an-ay  sizes  match  the  value  of  “n,”  the  variable  describing  the  number  of  tests  (line  17).  For  the  anay 
“censor,”  a  value  of  “zero”  is  given  to  indicate  a  test  nm  to  failure  whereas  a  value  of  “one”  is  given  to 
indicate  a  test  suspended  without  failure  (lines  24  to  29).  Array  “time”  (lines  31  to  45)  provides  the  test 
times.  The  subroutine  plvpercent  is  called  (line  56)  to  determine  the  relative  likelihood  ratio  for  a  speci¬ 
fied  10-percent  life.  By  calling  the  subroutine  repeatedly  within  a  do-loop  (lines  50  to  59),  data  for  a 
profile-likelihood  plot  are  calculated.  The  output  from  the  main  program  was  plotted  using  two  vertical 
scales  to  show  both  the  profile-lUcelihood  ratio  values  and  the  confidence  level  (fig.  3).  Table  III  lists 
some  profile-likelihood  ratios  and  the  con'esponding  confidence  level.  A  90-percent  confidence  interval 
for  the  10-percent  life  is  shown  graphically  on  figure  3.  The  lower  and  upper  bounds  can  be  determined 
more  precisely  by  running  the  program  again  with  an  appropriate  range  of  values  for  the  10-percent  life 
(figs.  4  and  5). 


TABLE  III.— CONFIDENCE  LEVELS 
AND  CORRESPONDING  PROFILE- 
LIKELIHOOD  RATIOS 


Confidence  level 

Ratio 

0.50 

0.7965 

.60 

.7018 

.70 

.5844 

.80  , 

.4399 

.85 

.3548 

.90 

.2585 

.95 

.1465 

.97 

.0811 

.99 

.0362 

Figure  3. — Profile-likelihood  ratios  and  confidence  levels  for  10-percent  life  for  data  of 
table  I. 
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Figure  4. — Profile-likellhood-based  confidence  interval  for  1 0-percent  life  for  data  of 
table  I  (lower  bound). 


Figure  5. — Profile-likelihood-based  confidence  interval  for  10-percent  life  for  data  of 
table  I  (upper  bound). 


The  subroutine  plvpercent  called  by  the  main  program  determines  the  shape  and  scale  factors  corre¬ 
sponding  to  the  largest  likelihood  value  by  an  iterative  search  method.  This  method  requires  that  maxi¬ 
mum  and  minimum  values  for  the  shape  parameter  be  provided  to  begin  the  search.  The  minimum  and 
maximum  values  used  in  this  example  were  hardcoded  in  the  subroutine  as  0.05  and  8.0  respectively. 
Users  who  wish  to  analyze  data  sets  that  require  inspection  of  shape  values  outside  the  just-stated  range 
need  to  modify  the  subroutine  (appendix  B). 


Source  Code  Listing  5 

Line  Source  Line  Microsoft  Fortran  PowerStation  Compiler.  Version  4.0 

1  I  a  main  program  to  calculate  likelihood  based  profile 

2  !  likelihood  curves  of  Weibull  distributions  from  censored  data 
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3  1 

4  !  this  program  presumes  that  a  percentile 

5  !  parameter  is  the  one  of  interest 

6  1  written  by  Tim  Krantz  On  2-13-2001 

7  implicit  none  !  declare  all  variables 

8  !  the  following  variables  represent  inputs 

9  1  they  are  hard-coded  in  this  main  program 

10  integer  n  '  the  number  of  tests 

11  real,  allocatable : :  time ( : )  I  the  test  times 

12  integer,  allocatable::  censor(:)  !  the  censoring  information 

13  integer  i  I  do  loop  counter 

14  real  blv, spblv, scblv 

15  real  pfrac 

16  real  tmin,tmax,dt 

17  n  =  15 

18  allocate (time (n) ) 

19  allocate (censor (n) ) 

20  1  HERE  WE  NEED  TO  SUPPLY  CENSORING  INFORMATION 

21  !  A  VALUE  OF  "0"  MEANS  THAT  A  FAILURE  OCCURRED 

22  !  A  VALUE  OF  "1"  MEANS  THAT  TEST  WAS  SUSPENDED  (CENSORED)  WITH  NO  FAILURE 

23  i  ENTER  THE  VALUES  AS  INTEGERS 

24  do  i=l,8 

25  censor(i)  =  0 

26  end  do 

27  do  i=9,n 

28  censor (i)  =  1 

29  end  do 

30  !  HERE  WE  NEED  TO  ENTER  THE  ARRAY  OF  TEST  TIMES 

31  time(l)=  40.2600 

32  time(2)=  53.9400 

33  time(3)=  59.8800 

34  time (4)=  67.6800 

35  time(5)=  95.7600 

36  time(6)=  134.220 

37  time(7)=  198.480 

38  time(8)=  256.200 

39  time(9)=  299.520 

40  time(10)=  301.560 

41  time(ll)=  303.660 

42  time(12)=  304.080 

43  time(13)=  305.800 

44  time(14)=  306.900 

45  time(15)=  335.400 

46  !  the  call  to  this  subroutine  checks  for  errors  in  the  censoring  array 

47  call  ccheck2 (n, censor) 

48  open (10, FILE=  "super_L10_profile.txt") 

49  pfrac  =  0.1 

50  tmin  =  4 . 

51  tmax  =  135. 

52  dt  =  (tmax-tmin)  /  100. 

53  tmin  =  tmin  -  dt 

54  do  i=l,101 

55  tmin  =  tmin  +  dt 

56  call  plvpercent (time, censor, n, tmin, pfrac, blv, spblv, scblv) 

57  write (10, 103)  tmin, spblv, scblv, blv 

58  103  format (4el6. 7) 

59  end  do 

60  stop 

61  end 


Data  Set  2 

The  main  program  that  determined  the  profile  likelihood  and  profile-likelihood-based  confidence 
intervals  for  the  10-percent  life  for  the  data  of  table  II  is  presented  below  as  source  code  listing  6.  The 
main  program  mirrors  that  of  source  code  5.  The  output  from  the  main  program  was  plotted  using 
two  vertical  scales  to  show  both  the  profile-likelihood  ratio  values  and  the  confidence  level  (fig.  6). 

A  90-percent  confidence  interval  for  the  10-percent  life  is  shown  graphically  on  figure  6. 
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Figure  6. — Profile-likelihood  ratios  and  confidence  levels  for  1 0-percent  life  for  data  of 
table  II. 


The  subroutine  plvpercent  called  by  the  main  program  determines  the  shape  and  scale  factors  corre¬ 
sponding  to  the  largest  likelihood  value  by  an  iterative  search  method.  This  method  requires  that  maxi¬ 
mum  and  minimum  values  for  the  shape  parameter  be  provided  to  begin  the  search.  The  minimum  and 
maximum  values  used  in  this  example  were  hardcoded  in  the  subroutine  as  0.05  and  8.0  respectively. 
Users  who  wish  to  analyze  data  sets  that  require  inspection  of  shape  values  outside  the  just-stated  range 
need  to  modify  the  subroutine  (appendix  B). 


Source  Code  Listing  6 

Line  Source  Line  Microsoft  Fortran  PowerStation  Compiler.  Version  4.0 

1  I  a  main  program  to  calculate  likelihood  based  profile 

2  1  likelihood  curves  of  Weibull  distributions 

3  I  from  censored  data 

4  I 

5  !  this  program  presumes  that  a  percentile  parameter  is  the  one  of  interest 

6  !  version  1 . 0 

7  !  written  by  Tim  Krantz  On  2-13-2001 

8  implicit  none  !  declare  all  variables 

9  !  the  following  variables  represent  inputs 

10  !  they  are  hard-coded  in  this  main  program 

11  integer  n  I  the  number  of  tests 

12  real,  allocatable : :  time  ( : )  !  the  test  times 

13  integer,  allocatable::  censor{:)  !  the  censoring  information 

14  integer  i  1  do  loop  counter 

15  real  blv, spblv, scblv 

16  real  pfrac 

17  real  tmin,tmax,dt 

18  n  ==  20 

19  allocate (time (n) ) 

20  allocate (censor (n) ) 

21  !  HERE  WE  NEED  TO  SUPPLY  CENSORING  INFORMATION 

22  !  A  VALUE  OF  "0"  MEANS  THAT  A  FAILURE  OCCURRED 

23  !  A  VALUE  OF  "1"  MEANS  THAT  TEST  WAS  SUSPENDED  (CENSORED)  WITH  NO  FAILURE 

24  do  i=l,17 

25  censor (i)  =  0 

26  end  do 

27  do  i=18,n 

28  censor (i)  =  1 

29  end  do 

30  1  HERE  WE  NEED  TO  ENTER  THE  ARRAY  OF  TEST  TIMES 
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31 

time  (1)  = 

7.08 

32 

time (2 ) = 

8.15 

33 

time (3 ) = 

15. 

52 

34 

time (4 ) = 

21 . 

34 

35 

time (5) = 

33. 

37 

36 

time (6 ) = 

35. 

70 

37 

time (7) = 

38. 

32 

38 

time (8 ) = 

40  . 

26 

39 

time (9) = 

40. 

74 

40 

time ( 10 ) = 

49. 

47 

41 

time (11) = 

51. 

70 

42 

time ( 12 ) = 

55. 

29 

43 

time (13 ) = 

59. 

17 

44 

time (14 ) = 

70. 

81 

45 

time (15) = 

127 

.07 

46 

time (16) = 

130 

.  95 

47 

time (17) = 

133 

.38 

48 

time (18) = 

300 

.00 

49 

time (19) = 

300 

.00 

50 

time (20) = 

300 

.00 

51  I  the  call  to  this  subroutine  checks  for  errors  in  the  censoring  array 

52  call  ccheck2 (n, censor) 

53  open  (10,  FILE=  "basel_L10_jprofile.txt") 

54  pfrac  =0.10 

55  write (6,*)  '  enter  min  value  ' 

56  read  (5,*)  tmin 

57  write (6,*)  '  enter  max  value  ' 

58  read (5,*)  tmax 

59  dt  =  (tmax-tmin)  /  100. 

60  tmin  =  tmin  -  dt 

61  do  i=l,101 

62  tmin  =  tmin  +  dt 

63  call  plvpercent (time, censor, n, tmin, pfrac, blv, spblv, scblv) 

64  write (10, 103)  tmin, spblv, scblv,blv 

65  103  format (4el6 . 7) 

66  end  do 

67  stop 

68  end 


Weibull  Plot  With  Confidence  Intervals 

Profile-likelihood-based  confidence  intervals  can  be  determined  for  any  percentile  of  interest  using 
main  programs  similar  to  source  code  listings  5  and  6.  The  variable  “pfrac”  determines  the  percentile  to 
be  investigated.  For  example,  using  “pfrac  =  0.30”  one  can  calculate  data  for  the  30-percent  life.  Collec¬ 
tions  of  such  confidence  intervals  were  calculated,  and  the  results  are  provided  in  figures  7  and  8.  These 
two  figures  are  plotted  using  Weibull  coordinates  so  that  the  Weibull  cumulative  distribution  frequency 
will  plot  as  a  straight  line.  The  test  data  points  are  plotted  at  the  positions  of  exact  median  ranks  (Jaquelin, 
1993). 
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Cycles  to  failure,  millions 


Figure  7. — Weibull  plot  for  data  of  table  I  with  profile-likelihood-based  confidence 
intervals.  Intervals  shown  are  pointwise  90-percent  confidence  intervals  of  percentile 
estimates. 


Cycles  to  failure,  millions 

Figure  8. — Weibull  plot  for  data  of  table  II  with  profile-likelihood-based  confidence 
intervals.  Intervals  shown  are  pointwise  90-percent  confidence  intervals  of  percentile 
estimates. 
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Summary  of  Results 


Softwaie  has  been  developed  to  analyze  gear  fatigue  test  data.  The  source  code  listings  discussed 
provided  examples  illustrating  how  to  define  the  test  data  to  be  analyzed  and  how  to  call  the  appropriate 
subroutines  to  analyze  the  data.  The  software  can  be  used  to  determine 

1.  Maximum  likelihood  estimates  of  the  Weibull  distribution 

2.  Data  for  contour  plots  of  relative  likelihood  for  two  parameters 

3.  Data  for  contour  plots  of  joint  confidence  regions 

4.  Data  for  the  profile  likelihood  of  the  Weibull  distribution  parameters 

5.  Data  for  the  profile  likelihood  of  any  percentile  of  the  distribution 

6.  Profile-likelihood-based  confidence  intervals  for  parameters  and/or  percentiles  of  the  distribution 

Contour  plots  of  relative  likelihood  (not  illustrated  in  this  document)  can  be  obtained  from  the  output 
files  as  provided  by  source  code  listings  3  and  4.  Profile-likelihood  calculations  for  Weibull  distribution 
parameters  were  not  illustrated  directly.  A  point  for  the  profile-likelihood  plot  of  the  shape  parameter  can 
be  determined  using  subroutine  wlratio  (appendix  B).  A  profile  likelihood  for  the  scale  parameter  can  be 
obtained  using  a  source  code  in  the  manner  of  listings  5  and  6  and  analyzing  for  the  63.2  percentile 
(“pfrac”  =  0.632). 
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Appendix  A 

Verification  of  Software 


The  software  was  verified  by  reproducing  results  published  by  Meeker  and  Escobar  (1998).  The  main 
programs  to  produce  the  figures  in  this  appendix  (figs.  9  to  14)  were  similai'  to  the  source  code  listings 
found  in  the  main  body  of  this  document. 


Figure  9. — Contours  of  equal  likelihood  ratios.  Weibull  shape  parameter  p  =  1/a 
and  Weibull  scale  parameter  nn  =  exp  ((jl).  In  the  manner  of  figure  8.1  from  Meeker 
and  Escobar  (1 998). 


9.8  10.0  10.2  10.4  10.6  10.8 


Figure  10. — Contours  of  equal  joint  confidence  levels.  Weibull  shape  parameter 
p  =  1/a  and  Weibull  scale  parameter  =  exp  (jjl).  In  the  manner  of  figure  8.4 
from  Meeker  and  Escobar  (1 998). 
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Confidence  level  Confidence  level 


Profile  likelihood 


6000  8000  10  000  12  000  14  000  16  000  18  000  20  000 

10-percent  life,  fo.i 

Figure  13. — Contours  of  equal  joint  confidence  levels  for  Weibull  shape  parameter 
and  10-percent  life,  fo.i  -  Weibull  shape  parameter  p  =  1/a.  In  the  manner  of 
figure  8.7  from  Meeker  and  Escobar  (1998). 


10-percent  life,  fo.i 

Figure  14. — Profile-likelihood  method  to  determine  confidence  interval  for  1 0-percent 
life,  fo.i  -  In  the  manner  of  figure  8.8  from  Meeker  and  Escobar  (1998). 
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Appendix  B 

Subroutine  Descriptions  and  User  Instructions 

Subroutine  CCHECK 

A  subroutine  to  do  “error”  checking  on  censoring  array,  all  values  should  be  either  1  or  0  integers. 
Usage 

subroutine  ccheck  (n,  censor) 

Arguments 

n  (integer)  size  of  aixay  “censor”  (input) 

censor  (integer  array)  array  of  dimension  [n],  censoring  information  (input) 

Comments 

If  all  elements  in  array  “censor”  equal  either  0  or  1 ,  no  action  is  taken. 

If  any  element  does  not  meet  the  criteria,  a  line  of  text  is  written  to  FORTRAN  units  numbered  6,  9, 
and  10,  and  then  execution  of  the  main  program  is  halted. 

Other  subroutines  called 

None 
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Subroutine  CCHECK2 


A  subroutine  to  do  enor  checking  on  censoring  array.  All  values  should  be  either  1  or  0  integers. 
Usage 

subroutine  ccheckl  (n,  censor) 

Arguments 

n  (integer)  size  of  an’ay  “censor”  (input) 

censor  (integer  airay)  array  of  dimension  [n],  censoring  information  (input) 

Comments 

If  all  elements  in  anay  “censor”  equal  either  0  or  1 ,  no  action  is  taken. 

If  any  element  does  not  meet  the  criteria,  a  line  of  text  is  written  to  FORTRAN  units  numbered  6  and 
10,  and  then  execution  of  the  main  program  is  halted. 

Other  subroutines  called 

None 
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Subroutine  FIXSHAPEMAXLL 


A  subroutine  to  calculate  the  maximum  likelihood  estimates  of  the  scale  parameter  for  a  two- 
parameter  Weibull  distribution,  for  a  presumed  shape  and  censored  data. 

Usage 

svbvoutimfixshapemaxll  (times,  censor,  n,  shape,  scale) 


Arguments 

(real  array) 
(integer) 
(integer  array) 
(real) 

(real) 


times 

n 

censor 

shape 

scale 


array  of  dimension  [n],  test  times  (input) 
size  of  an'ay  “censor”  (input) 

array  of  dimension  [n],  censoring  information  (input) 
presumed  Weibull  shape  parameter  (input) 
calculated  Weibull  scale  parameter  (output) 


Comments 

This  routine  is  based  on  the  code  published  by  Keats  and  Lawrence  (1997).  Although  arguments 
are  passed  as  single  precision,  calculations  within  the  subroutine  are  done  in  double  precision. 

For  array  “censor,”  a  value  of  0  indicates  a  test  to  failure  whereas  a  value  of  1  indicates  a  test  sus¬ 
pended  without  failure  (type  I  censoring). 

Other  subroutines  called 


None 
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Subroutine  PLVPERCENT 


A  subroutine  to  calculate  a  point  on  a  profile-likelihood  curve  for  a  Weibull  distribution  at  a  particu¬ 
lar  cumulative  distribution  function  (CDF)  percentile  of  interest  and  a  presumed  value  for  the  time  to 
failure  at  that  percentile.  This  routine  can  be  called  repeatedly  using  different  values  for  the  assumed 
time  to  failure  (atime)  to  generate  data  for  a  profile-likelihood  curve. 

Usage 

subroutine  plvpercent  (time,  censor,  n,  atime,  pfrac,  ratio,  spblv,  scblv) 

Arguments 

time 
censor 
n 

atime 
pfrac 

ratio 
spblv 
scblv 

Comments 

For  array  “censor,”  a  value  of  0  indicates  a  test  to  failure  whereas  a  value  of  1  indicates  a  test  sus¬ 
pended  without  failure  (type  I  censoring). 

Data  are  rescaled  within  the  subroutine  to  avoid  extreme  values. 

The  Weibull  shape  parameter  resulting  in  the  profile-likelihood  value  is  found  within  the  subroutine 
by  an  iterative  method.  The  shape  value  is  found  to  within  a  tolerance  of  0.01.  A  range  for  the  shape  must 
be  provided  to  start  the  iterative  process,  and  the  shape  parameter  resulting  in  the  profile-likelihood  value 
must  be  contained  within  that  range.  The  subroutine  uses  a  range  from  0.05  to  8.0.  These  values  could  be 
adjusted  to  fit  the  needs  of  a  particular  data  set,  but  one  might  encounter  numerical  problems  for  extreme 
scale  parameter  values  since  it  is  used  as  an  exponent  in  calculations.  The  routine  writes  a  line  of  text  to 
unit  number  6  if  the  routine  fails  to  converge  within  1000  iterations. 

Other  subroutines  called 

wmaxll 

wlratio 


(real  array) 
(integer  array) 
(integer) 

(real) 

(real) 

(real) 

(real) 

(real) 


art'ay  of  dimension  [n],  test  times  (input) 
array  of  dimension  [n],  censoring  information  (input) 
size  of  aiTay  “censor”  (input) 
assumed  time  to  failure  (input) 
assumed  percentile  of  CDF  corresponding  to 
value  of  atime  (input) 

profile  ratio  for  best  likelihood  value  [blv]  (output) 
shape  factor  corresponding  to  blv  (output) 
scale  factor  corresponding  to  blv  (output) 
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Subroutine  WBLINVERSE 


A  subroutine  to  calculate  an  inverse  Weibull  cumulative  distribution  function  (CDF)  value. 
Usage 

subroutine  wblinverse  (shape,  scale,  frac,  life) 

Arguments 


shape 

(real) 

Weibull  shape  parameter  (input) 

scale 

(real) 

Weibull  scale  parameter  (input) 

frac 

(real) 

Weibull  CDF  percentile  of  interest  (input) 

life 

(real) 

time  to  failure  for  provided  percentile  of  CDF 

Comments 

None 

Other  subroutines  called 

None 
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Subroutine  WFLSHAPE 


A  subroutine  to  calculate  Weibull  distribution  likelihood-based  profile  data  (likelihood  ratio). 

Usage 

subroutine  wjlshape  (z,  censor,  n,  shapemax,  shapemin) 

Arguments 

z 

censor 
n 

shapemax 
shapemin 

Comments 

For  aiTay  “censor,”  a  value  of  0  indicates  a  test  to  failure  whereas  a  value  of  1  indicates  a  test 
suspended  without  failure  (type  I  censoring). 

No  outputs  are  returned  to  the  calling  program.  However,  output  is  written  to  FORTRAN  unit 
number  9.  If  a  file  for  unit  number  9  was  not  yet  opened,  the  user  will  be  prompted  for  a  filename. 

The  likelihood  profile  ratio  is  calculated  for  41  values  of  the  shape  factor,  the  values  equally  spaced 
between  the  values  passed  as  shapemax  and  shapemin.  The  output  written  to  unit  number  9  consists  of 
two  columns  of  numbers  of  fixed  length,  where  the  numbers  are  the  shape  factor  and  the  calculated  likeli¬ 
hood  ratio. 

Data  are  rescaled  within  the  routine  to  avoid  extreme  values. 

Other  subroutines  called 

wmaxll 

fixshapemaxll 

wlratio 


(real  array)  array  of  dimension  [n],  test  times  (input) 

(integer  an'ay)  an*ay  of  dimension  [n],  censoring  information  (input) 

(integer)  size  of  array  “censor”  (input) 

(real)  maximum  shape  value  of  range  to  be  calculated  (input) 

(real)  minimum  shape  value  of  range  to  be  calculated  (input) 
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Subroutine  WJLONE 


A  subroutine  to  calculate  Weibull  distribution  likelihood-based  joint  confidence  ratios  and  cone- 
sponding  confidence  levels  to  be  plotted  as  contour  plots.  This  routine  presumes  that  the  Weibull  distribu¬ 
tion  is  parameterized  by  the  shape  parameter  and  10-percent  life. 

Usage 

subroutine  wjlone  (z,  censor,  n,  shapemax,  shapemin,  tenmax,  tenmin) 

Arguments 

z 

censor 

n 

shapemax 
shapemin 
tenmax 
tenmin 

Comments 

For  array  “censor,”  a  value  of  0  indicates  a  test  to  failure  whereas  a  value  of  1  indicates  a  test  sus¬ 
pended  without  failure  (type  I  censoring). 

No  outputs  are  returned  to  the  calling  program.  However,  output  is  written  to  FORTRAN  units 
numbered  9  and  10.  If  files  were  not  yet  opened,  the  user  will  be  prompted  for  filenames. 

The  likelihood  profile  ratios  and  confidence  levels  are  calculated  for  a  full  array  of  41x41  values 
(1681  lines)  of  the  shape  parameter  and  10-percent  life,  the  values  equally  spaced  between  the  ranges 
defined  by  the  values  passed  as  shapemax  to  shapemin  and  tenmax  to  tenmin,  respectively.  The  output 
written  to  unit  number  9  consists  of  three  columns  of  numbers  of  fixed  length,  the  numbers  being  the 
shape  parameter,  the  10-percent  life,  and  the  calculated  likelihood  ratio.  The  output  written  to  unit  number 
10  consists  of  three  columns  of  numbers  of  fixed  length,  where  the  numbers  are  the  shape  parameter,  the 
10-percent  life,  and  the  calculated  joint  confidence  level. 

Data  are  rescaled  within  the  routine  to  avoid  extreme  values. 

NOTE:  Confidence  levels  less  than  10  percent  are  written  to  the  output  file  as  equal  to  10  percent. 
Confidence  levels  greater  than  99.99  percent  are  written  to  the  output  file  as  equal  to  99.99  percent. 

Other  subroutines  called 

wmaxll 

wlratio 

c/i/iff  (IMSL  routine — see  IMSL  (1997);  IMSL  is  a  registered  trademark  of  Visual  Numerics,  Inc.) 


(real  array) 
(integer  anay) 
(integer) 

(real) 

(real) 

(real) 

(real) 


array  of  dimension  [n],  test  times  (input) 

array  of  dimension  [n],  censoring  information  (input) 

size  of  array  “censor”  (input) 

maximum  shape  value  of  range  to  be  calculated  (input) 
minimum  shape  value  of  range  to  be  calculated  (input) 
maximum  10-percent  life  to  be  calculated  (input) 
minimum  10-percent  life  to  be  calculated  (input) 
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Subroutine  WJLTWO 


A  subroutine  to  calculate  Weibull  distribution  likelihood-based  joint  confidence  ratios  and  corresponding 
confidence  levels  to  be  plotted  as  contour  plots.  This  routine  presumes  that  the  Weibull  distribution  is 
parameterized  by  the  shape  parameter  and  scale  parameter  (or  63.2-percent  life). 

Usage 

subroutine  wjltwo  (z,  censor,  n,  shapemax,  shapemin,  scalemax,  scalemin) 

Arguments 

z 

censor 
n 

shapemax 
shapemin 
scalemax 
scalemin 

Comments 

For  array  “censor,”  a  value  of  0  indicates  a  test  to  failure  whereas  a  value  of  1  indicates  a  test  sus¬ 
pended  without  failure  (type  I  censoring). 

No  outputs  are  returned  to  the  calling  program.  However,  output  is  written  to  FORTRAN  units 
numbered  9  and  10.  If  files  were  not  yet  opened,  the  user  will  be  prompted  for  filenames. 

The  likelihood  profile  ratios  and  confidence  levels  are  calculated  for  a  full  array  of  41x41  values 
(1681  lines)  of  the  shape  parameter  and  scale  parameter,  the  values  equally  spaced  between  the  ranges 
defined  by  the  values  passed  as  shapemax  to  shapemin  and  scalemax  to  scalemin,  respectively.  The  out¬ 
put  written  to  unit  number  9  consists  of  three  columns  of  numbers  of  fixed  length,  the  numbers  being  the 
shape  parameter,  the  scale  parameter,  and  the  calculated  likelihood  ratio.  The  output  written  to  unit  num¬ 
ber  10  consists  of  three  columns  of  numbers  of  fixed  length,  where  the  numbers  are  the  shape  parameter, 
the  scale  parameter,  and  the  calculated  joint  confidence  level. 

Data  are  rescaled  within  the  routine  to  avoid  extreme  values. 

NOTE:  Confidence  levels  less  than  10  percent  are  written  to  the  output  file  as  equal  to  10  percent. 
Confidence  levels  greater  than  99.99  percent  are  written  to  the  output  file  as  equal  to  99.99  percent. 

Other  subroutines  called 

wmaxll 

wlratio 

chidfilMSL  routine — see  IMSL  (1997);  IMSL  is  a  registered  trademark  of  Visual  Numerics,  Inc.) 


(real  airay) 
(integer  airay) 
(integer) 

(real) 

(real) 

(real) 

(real) 


array  of  dimension  [n],  test  times  (input) 

an  ay  of  dimension  [n],  censoring  information  (input) 

size  of  an'ay  “censor”  (input) 

maximum  shape  value  of  range  to  be  calculated  (input) 
minimum  shape  value  of  range  to  be  calculated  (input) 
maximum  scale  value  of  range  to  be  calculated  (input) 
minimum  scale  value  of  range  to  be  calculated  (input) 
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Subroutine  WLOGLV 


A  subroutine  to  calculate  the  log-likelihood  value  of  a  single  life  test.  Can  be  called  repeatedly,  with 
the  outputs  summed,  to  calculate  the  joint  log-likelihood  value  of  a  data  set. 

Usage 

subroutine  wloglv  (sc,  sh,  t,  fail,  llv) 


Arguments 


SC 

(real) 

Weibull  scale  parameter  (input) 

sh 

(real) 

Weibull  shape  parameter  (input) 

t 

(real) 

test  time  (input) 

fail 

(integer) 

flag  indicating  test  to  failure  of  test  censored  (input) 

llv 

(real) 

log-likelihood  value  (output) 

Comments 

For  variable  fail,  a  value  of  0  indicates  a  test  to  failure  whereas  a  value  of  1  indicates  a  test  suspended 
without  failure  (type  I  censoring). 

Other  subroutines  called 

None 
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Subroutine  WLRATIO 


A  subroutine  to  calculate  the  numerator  of  a  likelihood  ratio  of  a  data  set  using  presumed  Weibull 
shape  and  scale  parameters. 

Usage 

subroutine  wlmtio  (time,  censor,  n,  shape,  scale,  Iv) 

Arguments 

time 
censoi 
n 

shape 
scale 
Iv 

Comments 

For  an-ay  “censor,”  a  value  of  0  indicates  a  test  to  failure  whereas  a  value  of  1  indicates  a  test  sus¬ 
pended  without  failure  (type  I  censoring). 

Other  subroutines  called 

wloglv 


(real  airay)  array  of  dimension  [n],  test  times  (input) 

(integer  airay)  array  of  dimension  [n],  censoring  information  (input) 

(integer)  size  of  array  “censor”  (input) 

(real)  presumed  Weibull  shape  parameter  (input) 

(real)  presumed  Weibull  scale  parameter  (input) 

(real)  likelihood  value,  or  value  of  numerator  (output) 
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subroutine  WMAXLL 


A  subroutine  to  calculate  the  maximum  likelihood  estimates  of  the  scale  parameter  for  a  two- 
parameter  Weibull  distribution,  for  a  presumed  shape  and  censored  data. 

Usage 

subroutine  wmaxll  (times,  censor,  n,  shape,  scale) 

Arguments 

times 
n 

censor 
shape 
scale 

Comments 

This  routine  is  based  on  the  code  published  by  Keats,  Lawrence,  and  Wang  (1997).  Although  argu¬ 
ments  are  passed  as  single  precision,  calculations  within  the  subroutine  are  done  in  double  precision. 

For  array  “censor,”  a  value  of  0  indicates  a  test  to  failure  whereas  a  value  of  1  indicates  a  test  sus¬ 
pended  without  failure  (type  I  censoring). 

The  routine  finds  the  maximum  likelihood  value  for  the  shape  parameter  iteratively,  where  the  shape 
parameter  is  the  root  of  an  equation.  The  routine  is  considered  as  converged  if  the  step  size  of  the 
Newton-Raphson  method  is  less  than  0.0003.  If  the  routine  fails  to  converge  after  100  iterations,  a  line  is 
written  to  FORTRAN  unit  number  15. 

Other  subroutines  called 

None 


(real  array) 
(integer) 
(integer  an-ay) 
(real) 

(real) 


array  of  dimension  [n],  test  times  (input) 
size  of  array  “censor”  (input) 

array  of  dimension  [n],  censoring  information  (input) 
calculated  Weibull  shape  parameter  (output) 
calculated  Weibull  scale  parameter  (output) 
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